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1 Introduction 

One of the greatest challenges of modern astrophysics is understanding how 
galaxies, such as our Milky Way, form within the framework of the Big Bang 
cosmology. The current theory of structure formation, the extension of the 
Big Bang model called the Cold Dark Matter (CDM) scenario, predicts that 
galaxies form within extended massive dark matter halos built from smaller 
pieces that collided and merged, resulting in the hierarchy of galaxies, groups, 
and clusters observed today. The entire sequence of events is thought to be 
seeded by quantum fluctuations in the very early Universe and governed by 
mysterious "dark matter" which constitutes about 85% of all matter in the 
universe. Although the accurate properties of galaxies depend on complicated 
baryonic processes (radiative cooling, formation and evolution of stars, etc.) 
operating on small scales, we expect that overall spatial distribution of dark 
matter halos is closely related to the observed galaxy distribution. Here we 
present numerical simulations designed to study the formation, evolution and 
present day properties of such dark matter halos in different cosmological 
environments. 

In all simulations the spatially flat cold dark matter model with a cosmo- 
logical constant (/ICDM with = 0.3, Qa = 0.7, as = 0.9, and h = 0.7), 
favored by most current observations, has been assumed. 

2 Numerical simulations 

The Adaptive Refinement Tree (ART) A^-body code |8[ ^ was used to run 
all numerical simulation analyzed in this paper. In some of the simulations 
described below the code also included eulerian gasdynamics The code 
uses Adaptive Mesh Refinement technique to achieve high resolution in the 
regions of interests. The computational box is covered with a uniform grid 
which defines the lowest (zeroth) level of resolution. The code then reaches 
high force resolution by recursively refining all high density regions using 
an automated refinement algorithm. This creates an hierarchy of refinement 
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meshes of different resolutions, sizes, and geometries covering regions of inter- 
est. The refinement data structures and algorithms Q allow individual cubic 
cells to be refined. The shape of the refinement mesh can thus effectively 
match the geometry of the region of interest. This algorithm is well suited 
for simulations of a selected region within a large computational box, as in 
the simulations presented below. 

During the integration, spatial refinement is accompanied by temporal 
refinement. Namely, each level of refinement, I, is integrated with its own 
time step Aai = Z\ao/2', where Auq is the global time step of the zeroth 
refinement level. This variable time stepping is very important for accuracy 
of the results. As the force resolution increases, more steps are needed to 
integrate the trajectories accurately. In addition to spatial and temporal re- 
finement, simulations described below also use non-adaptive mass refinement 
to increase the mass (and correspondingly the force) resolution inside a spe- 
cific region. The multiple mass resolution is implemented in the following way 
. We first set up a realization of the initial spectrum of perturbations in 
such a way that initial conditions for a large number (1024'^) of particles can 
be generated in the simulation box. Initial coordinates and velocities of the 
particles are then calculated using all waves ranging from the fundamental 
mode k — 2tt /L to the Nyquist frequency k ~ 2tt / L x where L is the 

box size and N is the number of particles in the simulation. Particles out- 
side high-resolution regions are then merged into particles of larger mass and 
this process can be repeated for merged particles. The larger mass (merged) 
particle is assigned a velocity and displacement equal to the average velocity 
and displacement of the smaller-mass particles. Extensive tests of the code 
and comparisons with other numerical TV-body codes can be found in 0. 

With increasing number of particles and resolution (i.e., number of re- 
finement cells) the memory as well as the computing time requirement of 
the code increase. In practice, on nodes with 16 Gb of shared memory and 
up to 16 CPUs we are limited to simulations with < 256"^ particles (runtime 
is CPU-bound). The only way to overcome this problem is to use MPI to 
distribute computations accross nodes. We have developed two different MPI 
algorithms to handle larger simulations. In the first approach, we select spa- 
tially distinct objects of interest and simulate each of these objects on one 
node with high mass and force resolution whereas in the remaining part of 
the simulation box lower resolutions will be used. Each node uses then stan- 
dard ART with OpenMP on shared memory. The inter-node communication 
is minimal. In the second approach, one can divide the whole simulation box 
into a number of sub-boxes which will be handled by different nodes. Each 
node simulates its own sub-box with high mass and force resolution whereas 
other sub-boxes have lower mass and force resolution. Again one integration 
step is done by standard ART with OpenMP on each node. After each in- 
tegration step nodes exchange communications on positions, velocities and 
masses of particles that cross sub-box boundaries. 
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2.1 Constrained simulations of the local universe 

Here our goal is to perform simulations that match the observed local universe 
as well as possible. Namely, we are interested in reproducing the observed 
structures: the Virgo cluster, the Local Supercluster and the Local Group, in 
the approximately correct locations and embedded within the observed large- 
scale configuration dominated by the Great Attractor and Perseus-Pisces 
superclusters. 

An efficient algorithm for reconstructing the density and velocity fields 
from sparse and noisy data of redshift and velocity surveys is provided by 
the Wiener filter formalism and constrained realizations of gaussian random 
fields Here we use the MARK HI catalog The sample consists of 
~ 3400 galaxies and provides radial velocities and inferred distances with 
fractional errors ~ 17 — 21%. A detailed analysis of the large scale structure 
in the Local Universe reconstructed from the MARK III survey was presented 
in 111. 

Several simulations with increasing force and mass resolution in the re- 
gion around the Virgo Cluster were performed |l^ . The initial conditions 
for these simulations were set using multiple mass resolution technique. Us- 
ing z = output of a low-resolution run, we selected all particles within a 
sphere of 25/i~^Mpc radius centered on the Virgo cluster. The mass resolu- 
tion in the Lagrangian region occupied by the selected particles was increased 
and additional small-scale waves from the initial ylCDM power spectrum of 
perturbations were added appropriately |^ . For the two high-resolution sim- 
ulations, the particle mass in the Local Supercluster region is 8 and 64 times 
smaller than in the low-resolution simulation. The highest resolution simula- 
tion has a particle mass of 2.5 x 10^h~^MQ and the maximum formal force 
resolution was 2.4/i^^kpc in the Local Supercluster region. The results of 
both high-resolution simulations agree well with each other at all resolved 
scales. 

All major structures (the Local Supercluster, Great Attractor, Perseus- 
Pisces supercluster, and Coma cluster) observed within 100/i^^Mpc around 
the Milky Way exist in the simulations. The positions and morphology of 
these structures is, of course, fairly well dictated by the constraints imposed 
on the initial conditions. The Local Supercluster is an elongated structure 
which extends over ^ 40/i^^Mpc along the SGX axis. There is a low-density 
"bridge" (of overdensity just above the average density), which connects the 
Local Supercluster with the Perseus-Pisces Super-cluster. There is also an 
even weaker filament connecting the Local Super-cluster with the Great At- 
tractor. 

Just as in the real Universe, the Local Group is located in a weak fil- 
ament extending between the Virgo and Fornax clusters. This filament is 
a counterpart of the Coma-Sculptor "cloud" in the distribution of nearby 
galaxies. Figure ^ shows a zoom-in view of the immediate environment of 
the simulated Local Group. Note that the structures at these scales are only 
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Fig. 1. Density (contours corresponding to overdensities of 1, 2, and 3) and velocity 
(arrows) fields smoothed with a Gaussian filter of 0.7/i~^Mpc smoothing length 
around the Local Group. The slice shown has a size and thickness of 15/i^^Mpc 
and 5/i^^Mpc, respectively, and is centered on the supergalactic plane {SGZ — 0). 
The length of the thick arrow in the top right corner corresponds to a velocity of 
500 km/s. The velocities are plotted in the Virgo cluster rest frame. 

weakly affected by constraints imposed on the initial conditions. Several pos- 
sible counterparts to existing objects (e.g., the MW and M31, M51, NGC253) 
are marked, but their existence is largely fortuitous. As can be seen in this 
figure, the simulated Local Group is located in a rather weak filament extend- 
ing to the Virgo cluster. This filament borders an underdense region visible 
in the right lower corner of Figure |l|, which corresponds to the Local Void 
in the observed distribution of nearby galaxies. We are now carrying out a 
series of high-resolution simulations of such voids (see below) . The velocity 
field around the Local Group is rather quiet, in good agreement with ob- 
servations. The peculiar velocity field in the Local Void exhibits a uniform 
expansion of matter out of this underdense region, while velocities between 
the Local Group and the Virgo (upper half of Fig. |l|) show a coherent flow 
onto the main body of the Local Supcrcluster. 

2.2 Galaxy clusters 

Within a low mass resolution simulation (128^ particles, rUpart = 2.0 x 
10^°/i~"'^Mq) we have identified 15 clusters of galaxies with masses above 
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Fig. 2. Merger event between two 6 x 10 M0 clusters. Particles belonging to one 
of the two clusters before merging are plotted with red and blue dots, respectively. 
Yellow lines correspond to density contours. 

1.0 X IO^^H-^Mq. From this set we have selected 8 candidates with differ- 
ent masses and merging histories and added 5 smaller clusters/groups for 
load balance. We then re-simulated the clusters with higher mass resolution. 
With particle masses of 3.2 x 10^ Mq a typical cluster and its environ- 
ment contains more than 1 million particles. The formal force resolution with 
9 refinement levels was l/i~^kpc. Halos with masses above 3.0 x IO^'^H'^Mq 
are well resolved. A typical cluster contains more than 150 such halos. The 
simulations were done using an MPI version of the ART code where each of 
8 nodes followed the evolution of one or two clusters. We are now carrying 
out simulations of these clusters including gasdynamics. 

The dynamical evolution of a typical major merger between galaxy clus- 
ters can be seen in Fig. ^. The figure shows three snapshots of a cluster 
taken at z = 0.8, z = 0.6, and z — 0.2. Although the cluster does not show 
significant substructure in its density contours during the merger, it is far 
from a relaxed state. Most of the energy transfer takes place during the first 
encounter (between z = 0.8 and z = 0.6), but the merger remnant does not 
reach virial equilibrium until it loses all the information about the initial con- 
ditions (i.e. the particles are well mixed). For this event, virialisation occurs 
at 2; ~ 0.2 (~ 4 Gyr after the first encounter). During the first stages of clus- 
ter formation, the characteristic time between major mergers can be shorter 
than the relaxation time. High-redshift clusters may thus be in general far 
from the virial equilibrium inside their formal virial radius. 

2.3 Voids 

Cosmological simulations predict many more small DM halos than the ob- 
served number of satellites around the Milky Way and Andromeda galaxies. 
Do we have the same problem for dwarf galaxies in voids? One naively ex- 
pects a large number because the Press-Schechter mass function steeply rises 
with declining mass. In contrast, it seems that observations are failing to find 
a substantial number of dwarf galaxies inside voids (e.g. ||Tl|, |3|). However, the 
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Fig. 3. The inner region of a spherical void at 2: = 0. 

situation is complicated because it is very difficult to detect dwarf galaxies, 
many of them are expected to have low surface brightness. 

In order to study the formation of large voids, the simulation box should 
be sufficiently large; we use a cube of size 80/i~^Mpc. On the other hand, we 
are interested in the formation of small structure elements inside voids, for 
which we need highest possible mass resolution. Therefore, we use particles 
with large masses to follow the evolution of large scale structures and particles 
with small masses to follow the evolution of structure within one spherical 
void. 

Within a low mass resolution simulation (128^ particles, nipart = 2.0 x 
10^^h~^MQ) we have identified 8387 galactic halos with masses > 2.0 x 
10^^h~^M(7). This corresponds to a mean distance of about 4/i~^Mpc between 
halos. We then searched voids in the distribution of these galactic halos by 
constructing the minimal spanning tree using halo positions. The minimal 
spanning tree was then used to search for the point in the simulation box 
which has the largest distance ri to the set of halos. We identified this point 
as the center of spherical void with the radius of ri. Excluding that void we 
were searching again for the point with the largest distance to the set and 
thus found the second largest void and so on. The algorithm is similar to that 
used by §. 

With the algorithm described above we find spherical voids in the halo 
distribution which do not contain any halo with a mass larger than 2.0 x 
lQ^^h~^MQ; such halos by definition lie on the border of the void. After 
finding voids in the low-resolution simulation, we re-run the simulation with 
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much higher mass resolution inside the voids (mpart = 4.0 x IO'^H^^Mq), 
which allows identification of objects with masses larger than 10^/i~^Mq. 
By construction, the voids do not have halos with masses larger than Mb = 
2.0 X 10^^h~^MQ. Five voids have been identified and resimulated, their radii 
are r^oid = 11.6, 10.8, 9.4, 9.1, 9.1 ft,~^Mpc. Inside the voids the matter 
density is typically a factor of 10 smaller than the mean density, but void 
properties exhibit large differences. Some voids are very isolated: the density 
within a sphere with radius 30 h~^Mpc centered on one of the void centers 
reaches only 2/3 of the mean density. On the other hand, the most prominent 
structure of the simulation, a galaxy cluster of 2 x 10^^h~^Mpc, is bordering 
one of the other voids. The density outside of this void rapidly increases. The 
mass function in voids is about an order of magnitude lower and its shape is 
different than that of the field galaxies |^ . 

In Fig. H we show the inner 20 ft,~^Mpc of a spherical void of radius 21.6 
h~^Mpc. In this void we found more than 50 halos with circular velocity 
Vc > 50 km/s and more than 600 halos with Vc > 20 km/s. There is a certain 
spatial mass segregation among halos. Typically, more massive halos tend to 
be situated in the outer part of the void. The largest halos in the plot are 
actually dwarf-size halos with circular velocities of ~50km/s. The void is far 
from being empty and boring: it has a complex structure with numerous long 
filaments and small sub-voids. Visually it resembles the large-scale structure 
of the Universe, but everything in this plot is hundreds and thousands times 
less massive than in "normal" configurations of interconnected superclusters 
and filaments. 



2.4 Hydro simulations 

We have started to repeat |jl^ many of the simulations described above in- 
cluding eulerian shock-capturing gasdynamics and physical processes such as 
radiative cooling and stellar feedback. Inclusion of gasdynamics, although sig- 
nificantly increasing computational and memory demands of simulations, will 
allow us to address a much wider range of questions and more robust compar- 
ison with observations. For example, simulations of the Local Supercluster 
regions with cooling and starformation will allow us to study distribution 
and properties of galaxies in the simulations as a function of their luminosity 
and color. High-resolution simulations of the local voids will allow studies 
of spatial distribution of Ly a absorbers and their connection to galaxies. 
Simulations of nearby galaxy clusters (e.g., Virgo, Fornax, Coma) formed 
in realistics large-scale environments will allow for detailed object-to-object 
comparison of properties and should give us good insights into physical pro- 
cesses operating within intracluster medium. 

As an example, we present here simulations of the Virgo cluster in the 



context of constrained simulations described in § 2.1. The lagrangian region 



within five virial radii (at ^ = 0) around the cluster was resimulated with 
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Fig. 4. The projected X-ray intensity and emission-weighted temperature maps of 
the Local supercluster in galactic coordinates. The Virgo cluster is the brightest 
object near the center (which corresponds to the North Galactic Pole). 

mass resolution 8 times higher (particle mass of 3.1 x IO^/i^^A/q) than that 
of the Local Supercluster simulation and spatial resolution of ~ l/i~^kpc in 
the central regions of the cluster. Figure ^ shows the sky projection of the 
X-ray intensity and emission- weighted temperature of gas in the simulations 
in galactic coordinates (the center of the polar projection corresponds to the 
North Galactic Pole, NGP). The Virgo cluster is the brightest object located 
near the NGP, very close to the actual location of the real Virgo cluster. Fig- 
ure ^ shows radial profiles of dark matter and gas density, as well as profiles 
of gas temperature and entropy for two simulations: one that included only 
adiabatic gasdynamics, and the other in which gas was preheated with the en- 
ergy of 1.5 keV per gas particle at z = 3. The preheating very roughly models 
the possible effect of galactic winds at high-redshifts. The figure shows that 
preheating increases entropy of the gas thereby lowering its density in the 
cluster core and increasing the overall gas temperature. It lowers the X-ray 
luminosity of the cluster by a factor of 8, which brings it in good agreement 
with the observed luminosity-temperature relation. Note, however, that it 
does not change the shape of the temperature profile. The proximity of the 
Virgo cluster makes it one of the best spatially resolved clusters, with the 
temperature mesurements well outside the virial radius. Figure ^ compares 
the observed temperature measurements around the center of the Virgo clus- 
ter |]l^ as a function of projected radius to the corresponding measurements 
in simulations (adiabatic and preheating). We constructed projected map of 
emission-weighted temperature of the simulated Virgo cluster using 4' x 4' 
pixels, each pixel represented by a point in the figure. The figure shows that 
neither simulation can match the much shallower temperature distribution in 
the central regions of the Virgo cluster. This indicates that preheating alone 
cannot be the whole story and other processes, such as cooling and/or central 



Simulations of the Local Universe 



9 




» Adiabatic 
o Preheating: 

1,5 keV / par. @ z = 3 



m\ I I I llll| 1 



0.01 0.1 1 

r [h-' Mpc] 



^ 2.5 

M b O o O ° 

2 



0.01 0.1 1 10 

r [h"' Mpc] 



Fig. 5. Radial profiles of dark matter and gas density, gas temperature, and gas 
entropy of the simulated Virgo cluster. The dotted line in the DM density panel 
shows the NFW profile with the concentration of 7. 

heating by AGNs can be important. We will be exploring the effects of these 
processes in future simulations. 



3 Summary 

The adaptive refinement tree code is a useful tool to study cosmological struc- 
ture formation on different scales and with different resolutions. It runs well 
on a variety of platforms with shared and distributed memory. The simula- 
tions described here have been performed on the Hitachi of LRZ Munich, the 
small development Hitachi at the AIP, the IBM SP of the Potsdam Institute 
for Climate Impact Research, and the IBM SP of NERSC Berkeley. 
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